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Many neuronal systems and models display a certain class of mixed mode oscillations (MMOs) consisting 
of periods of small amplitude oscillations interspersed with spikes. Various models with different underlying 
mechanisms have been proposed to generate this type of behavior. Stochastic versions of these models can 
produce similarly looking time series, often with noise-driven mechanisms different from those of the deter- 
ministic models. We present a suite of measures which, when applied to the time series, serves to distinguish 
models and classify routes to producing MMOs, such as noise-induced oscillations or delay bifurcation. By 
focusing on the subthreshold oscillations, we analyze the interspike interval density, trends in the amplitude 
and a coherence measure. We develop these measures on a biophysical model for stellate cells and a phe- 
nomenological FitzHugh-Nagumo-type model and apply them on related models. The analysis highlights the 
influence of model parameters and reset and return mechanisms in the context of a novel approach using noise 
level to distinguish model types and MMO mechanisms. Ultimately, we indicate how the suite of measures 
can be applied to experimental time series to reveal the underlying dynamical structure, while exploiting 
either the intrinsic noise of the system or tunable extrinsic noise. 

PACS numbers: 87.18.Tt Noise in biological systems; 05. 45. -a Nonlinear dynamics and chaos; 05.45.Tp Time 
series analysis 



I. INTRODUCTION 

Mixed mode oscillations (MMOs), composed of well- 
defined periods of subthreshold oscillations (STOs) sepa- 
rating spikes, appear in biological rhythms, neural dy- 
namics, chemical oscillations, and network oscillators, 
with many models proposed to capture this phenomenon 
in these and other applications^. One question in ap- 
plications is the significance of the STOs - that is, do 
the interspike intervals (ISIs) and STOs encode impor- 
tant information for the system (e.g., Refs. Isl or[23l). An 
element in the answer is identifying the mechanisms that 
drive the STOs. A challenge in this identification is that 
noisy time series generated by structurally different mod- 
els can appear hardly distinguishable from each other, 
even quantitatively. Experimentally observed behavior 
can be captured with different types of models with dif- 
ferent noise levels, so that model identification and cali- 
bration requires consideration of several classes of models 
over wide parameter ranges. 

Previous studies indicate that a minimum of three de- 
grees of freedom is necessary to produce MMOs in deter- 
ministic models (see examples in Ref. , while stochastic 
van der Pol and FitzHugh-Nagumo-type models illustrate 
how noise can drive MMOs in simple 2D models^i^l. 
Similarities in MMO signals often stem from an under- 
lying Hopf bifurcation and canard structure in the sys- 
tem^ i^^'^° . In periods where the system is in some sense 
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close to the Hopf bifurcation, STOs can be both sup- 
ported and disrupted by noise. In these models it is pos- 
sible to choose parameters so that the MMOs are driven 
by the deterministic dynamics, while for other parame- 
ter choices the appearance of MMOs relies on stochastic 
fluctuations. 

While on the one hand noise can make it difficult to dis- 
tinguish between different models with various routes for 
producing stochastic MMOs, it is nevertheless possible 
to use noise as a way to identify underlying mechanisms. 
In this paper we compare a suite of measures at different 
noise levels to extract the structure of the model, and 
thus the mechanism for the STOs and MMOs. 

The choice of measures is partly motivated by the pres- 
ence of multiple time scales, an important feature that 
allows both well-defined periods of STOs and more than 
one mechanism for MMO generation. As described in 
more detail in the next section, the common substruc- 
ture of the models we consider corresponds to a subsys- 
tem with a slowly varying control parameter. Depending 
on the combination of the other model parameters, that 
control parameter can slowly vary through a Hopf bifur- 
cation, in which case the STO-spike transition exhibits 
the dynamics of bifurcation delay (Ref. 6 and references 
therein). For other parameter combinations correspond- 
ing to quiescence in the deterministic system, noise ex- 
cites subthreshold coherent oscillations with a frequency 
near that of the Hopf bifurcation. The weak damping of 
these oscillations is due to the presence of multiple time 
scales, observed for parameters near a Hopf bifurcation. 
These noise-induced STOs are just one type of oscillation 
appearing in the context of coherence resonance (CR). In 
its broadest sense CR is observed when a range of noise 
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levels excites coherent oscillations in a system that is qui- 
escent without noise, with a maximum coherence at an 
optimal noise level. This phenomenon is most commonly 
cited for transitions between steady equilibria with large 
excursions, such as relaxation oscillations^^'^^, but in this 
paper it refers to the noise- induced STOs near a HB, also 
exhibiting CR-type behavior. 

Whether the STOs appear via slow passage through a 
Hopf bifurcation (SPHB) or via a CR-type mechanism, 
the impact of the noise on the MMOs is concentrated in 
the interspike interval (ISI) where the STOs arise and 
eventually transition to a spike. It is no accident that 
this is also the interval where the slow dynamics are 
most prominent, and it is well known that noise can 
have a significant impact in such intervals^^i^. The 
suite of measures used in this paper focuses on the STOs 
and the ISI behavior, capitalizing on previous analytical 
and computational studies of noise sensitivity of SPHB'' 
and CR-driven oscillations^7'36^ These previous analy- 
ses show how noise can enhance or suppress the STOs 
through interactions with the fast-slow dynamics inher- 
ent in the MMOs, and these characteristics are captured 
by the behavior of the suite of measures. Focused on the 
noise sensitive STOs in the ISI, these measures can iso- 
late differences between models and identify mechanisms 
for the MMOs as noise levels and bifurcation and control 
parameters vary. 

An additional advantage of these measures is that they 
are based on features that are easy to analyze in time 
series data, making them amenable for use on experi- 
mental and simulated data. Furthermore, recognition of 
the multiple time scales of the MMOs suggests ways to 
experimentally introduce fluctuations into the system as 
a tunable extrinsic noise. When scaled appropriately by 
considering the ISI dynamics, this extrinsic noise can be 
used to mimic the effects of other intrinsic noise sources. 
Then this introduction of noise, combined with the suite 
of measures, can be applied both to identify the likely 
mechanism for MMOs and limit the variety of models 
one must consider for calibration. 

Summary of results 

We use a combination of amplitude trend, ISI den- 
sity, and coherence measure based on power spectral den- 
sity (PSD) to classify and differentiate mechanisms for 
MMOs. We compare phenomcnological and physiological 
models to explore the types of STO and MMO behaviors 
that can be captured by simplified or reduced stochastic 
models over a range of noise levels. We then analyze these 
results within intra-model comparisons to get a deeper 
understanding of those MMO characteristics related to 
model structure, those dictated by deterministic dynam- 
ics, and those driven by noise. We focus particularly on 
cases where CR and/or SPHB are key mechanisms. In 
the conclusion we give a detailed list of identifying char- 
acteristics of MMO mechanisms obtained from this suite 
of measures. Understanding the dynamical features re- 
veals the sensitivity in model calibration on bifurcation 
structure and on model choice, and illustrates ways that 



noise can be used to help to differentiate between models. 

For STOs of the CR-type we observe longer tails in 
the ISI density together with stronger peaks in the co- 
herence measure vs. noise, and no trend in the aver- 
age STO amplitude leading up to the spike. In contrast, 
SPHB-dominated STOs have different behaviors in all 
three measures. 

Simplified low dimensional integrate and fire (IF) mod- 
els can be used to capture some STO and ISI behavior 
over a range of noise levels through adjustment of the 
speed of the control variable and reset. However, they 
can not mimic differences observed from underlying bi- 
furcation structure, particularly in the case where a com- 
bination of CR and SPHB are at play. 

Models with voltage-dependent control variables can 
have a rich variety of underlying deterministic behav- 
ior, allowing possibilities for both CR- and SPHB-driven 
MMOs. These differences can drive significantly differ- 
ent stochastic behavior in the ISI and PSD behavior as 
compared to simple IF models. 

Related to this last observation, our analysis also il- 
lustrates how refractory dynamics, dynamics following 
the spike at the initial stage of the ISI, can influence ISI 
density for increasing noise, and disrupt the coherence 
of the STOs in MMOs. The influence of the reset or 
return mechanism, as well as the underlying determinis- 
tic behavior, observed in the measures considered here is 
consistent with the influence of these modeling aspects 
observed for noise-driven clusters of spikes, that is, re- 
peated spikes without STOs in the ISI^^. 

The article is organized as follows: In Section |TT] we 
first give a brief overview of two mathematical models in 
the literature - one physiological, the other phenomcno- 
logical - that have been used to explain MMOs in neural 
systems, and discuss different mechanisms for MMOs in 
these stochastic systems. In the same section we present 
a number of different computational measures for analyz- 
ing important characteristics of the time series in order to 
identify the underlying model structure. In Section lllll we 
generate and classify MMOs with the physiological model 
and find matching MMOs generated by the phenomcno- 
logical model. We then apply the measures to these time 
series and provide a set of features that - in combina- 
tion - make the time series distinguishable. Sees. IIVI 
and |V] focus on the comparison of time series generated 
by the same model with varying parameters, to further 
highlight differences in the computational measures that 
appear with different model or bifurcation structure. We 
also highlight the dramatic effect of weak noise on differ- 
ent families of MMOs that are close in parameter space. 
In Sec. IVII we outline how tunable extrinsic noise can be 
introduced to mimic the effect of intrinsic noise, based 
on the inherent multiple time scales and the fact that 
the noise impact is concentrated in the ISI. We also illus- 
trate the applicability of our analysis to another MMO- 
generating model. 
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II. MODELS AND MEASURES 

We review the setup and dynamics of two models 
whose STO dynamics will be analyzed in this paper: A 
detailed biophysical model for MMOs in stellate cells 
('SC model) as well as a modified FitzHugh-Nagumo 
model ('MFN'). The two main MMO-generating mecha- 
nisms in these noisy models are slow passage through a 
Hopf bifurcation (SPHB) and coherence resonance (CR). 

A. Biophysical model - SC 

We consider a biophysical model for the MMOs ob- 
served experimentally in the layer II stellate cells of the 
medial entorhinal cortex introduced in Ref. [l]. It con- 
sists of a voltage equation that includes the three usual 
Hodgkin-Huxley currents for sodium, potassium and a 



leak current together with a persistent sodium current 
and a hyperpolarisation-activated current. Together with 
the dynamical equations for six gating variables, this 
model is a system of seven coupled ordinary differential 
equations ('7DSC'), given in App. [Bl 

Since we focus on the STO dynamics, for most of this 
paper we analyze a reduced, three-dimensional version 
(3DSC) of the full model, introduced in Ref. HO, consist- 
ing of only the equations for the voltage and two gating 
variables rf and (In the equation for r^, we use the 
term from the original work^ (see footnote in Ref. ,3(|).). 
This 3DSC model provides a good approximation of the 
ISI dynamics of the full 7D model?-. Throughout this 
article, we refrain from giving units in the equations or 
for parameters (see App. [B] for the original units and pa- 
rameter values). The three equations considered for the 
reduced model are: 



lapp — Gl{V — El) — Gp 



'exp(-^) 



- 0.15V2Dr]it) I (V - £;Na) - G,.(0.65r/ + 0.35r,){V - Eh) 



ff = [l/[l + exp((y + 79.2)/9.78)] -r/]/[0.51/[exp((F- 1.7)/10) + cxp(-(T/ + 340)/52)] + 1] 
j-s = + exp {{V + 71.3)/7.9)] - r,] / [5.6/[exp {{V - 1.7)/14) + cxp {-{V + 260)/43)] + 1] . 



(1) 
(2) 
(3) 



Eq. [T] contains a noise term 'q{t) representing white 
noise, corresponding to fluctuations in the persistent 
sodium current. According to Ref. [H, this current gives 
the main stochastic contribution in stellate cells, due to 
the low number of channels. The model (Eqs. [lH3|) is 
treated as an integrate and fire (IF) model, where spik- 
ing and refractory dynamics of the full (7D) model are 
replaced by a reset value once the trajectory of the system 
crosses a threshold defined in terms of V . The value of 
reset is chosen such that the trajectory is placed close 
to a post-spike value of the full deterministic system 
(F = — 80, r/ = = as in Ref. [13) • The deterministic 
[D — 0) version of this reduced model has been analyzed 
in detail^i^i^, and there have been also studies of the 
noisy (D ^ 0) casei^i^. 

Both the full 7D and the 3DSC deterministic sys- 
tem produce attracting states of either a steady state, 
MMOs, or tonic spiking, depending on /app- In the 3DSC 
model these different long time patterns are observed for 
/app < /app.H ~ -2.575, -2.575 < /app < -2.241, and 
/app ^ —2.241, respectively. These values are slightly 
shifted in the 7D version of the model, e.g., the Hopf bi- 
furcation value /app.H separating steady state solutions 
from MMOs lies at /Jpp^ « -2.702 when using the al- 
ternative Eq. IB9I for the dynamics of Ts (see App. IB|). 

As discussed in Ref. |30|, the 3DSC system can be 
viewed as a 2D subsystem in V-rj with slowly varying 



control parameter . In the nondimensionalized version 
of the 3DSC model analyzed in Ref. Ill], the ratio of time 
scales for V and rf is approximately 0.02, with the time 
scale for even slower than that of by a factor of 0.3. 

From that viewpoint. Fig. [T] shows the relevant part 
of the bifurcation diagram (with fixed /app) of the un- 
derlying 2D subsystem with treated as the bifurcation 
parameter. For the parameters considered, there is a sub- 
critical Hopf bifurcation at r^^H and a canard transition 
at rs,c {tsM ~ 0.08437, r^.c « 0.08271 for /app = -2.45 
and r,,H « 0.09241, r,^c ~ 0.09078 for /app = -2.58 - 
obtained with XPFAUT^O). 

In the stochastic 3DSC model {D > 0), MMOs are 
observed for a broader range of /app values than listed 
above, generated by two mechanisms: a slow passage 
through a Hopf bifurcation and coherence resonance. For 
intermediate values of /app and noise, mixtures of both of 
these mechanisms are observed. One of the purposes of 
this paper is to provide measures that characterize and 
identify these different mechanisms in various settings. 
The top panel of Fig. [T] shows the stochastic dynam- 
ics of the 3DSC model when the deterministic dynam- 
ics corresponds to a SPHB, or delay bifurcation. The 
dependent variable V is attracted to a steady state for 
Ts < Ts,¥i- As the variable slowly ramps through r^.u, 
the full system does not immediately make the transi- 
tion to spiking. Instead, there is a delayed transition. 
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as STOs increase gradually about the unstable steady 
state. It is well known from previous studies that for 
delay bifurcations, or more generally for certain systems 
with alternating slow/fast dynamics, certain character- 
istics of the slow dynamics are exponentially sensitive 
to noise. That is, for noise levels above those exponen- 
tially small in terms of the parameter of the slow time 
scale, the behavior of the slow dynamics is qualitatively 
changed by the noise, as has been discussed for exam- 
ples of delay bifurcations and resonance^"— and in 
other slow systemsi2ii2ii2i. In the example in Fig. [T] the 
stochastic model typically has shorter intervals of STOs 
as compared with the deterministic model, as even small 
noise typically drives faster transitions to spiking. 

The example in the bottom panel of Fig. [T] corre- 
sponds to STOs of the CR-type. The underlying de- 
terministic system has a steady state corresponding to 
"Ts = i^s,o < 1^3,11, so STOs are damped. In the stochas- 
tic system coherence resonance (CR) occurs when noise 
amplifies the STOs in the presence of weak damping (see 
Refs. [nl, [13, HO, m and references therein; note that 
this type of CR is different to what is described, e.g., 
in Refs. [l^ or [28h. These STOs have a frequency corre- 
sponding to that of the Hopf bifurcation at rs,H of the 
2D subsystem. Results in Ref. [H indicate that CR drives 
STOs in the range of noise levels where coherence is op- 
timal, that is, where the PSD is relatively narrow with 
significant power. The analysis in Ref. 3d shows that 
the amplification factor of the STOs is inversely propor- 
tional to the proximity to the Hopf bifurcation, so that 
the STOs are a prominent and prolonged feature in the 
ISI for Vs near r^^. Variability in the amplitude of these 
STOs yields a nontrivial probability of spiking, leading 
to frequent appearance of MMOs, even for /app < /app,H 
if noise is strong enough^. 

B. Modified FitzHugh-Nagumo model - MFN 

In Ref. a FitzHugh-Nagumo-type (FN) model is 
given by the following pair of coupled nonlinear stochas- 
tic differential equations: 

eu — u{u — a)(l — u) — v, (4) 

v = g{u-b) + V2D^{t). (5) 

We follow Ref. H and use e = 0.005, a = 0.9 and 
g{x) = 7x^ + 0.08(1 — e~^/° '^^), chosen to give a ratio 
of the time scales of the STOs and the spikes (relaxation 
oscillations) that is close to that typically observed exper- 
imentally and in the SC model. The parameter 6 is a con- 
trol or bifurcation parameter. We refer to this model as 
the modified FN model (MFN). Taking a ^ -I, D ^ 0, 
adding a constant term in Eq. 0] and a w-dependent term 
in Eq. [5] yields the general FN modeli^. With a = -1, 
g{x) = X, b = D = and a simple rescaling Eqs. |4]and[5] 
are the original van der Pol (vdP) model"'^'^. 

Fig. [5] shows the corresponding bifurcation diagram: a 
supercritical Hopf bifurcation with a canard transition to 
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FIG. 1. Examples of the dynamics of the 3DSC model (for 
units see App.[B]): trajectory in the underlying 2D bifurcation 
diagram and time series. Top: slow passage through a HB 
(see Class 1 later; 7app = -2.45, D = 10"^); bottom: STOs 
of the CR-type (see Class 2 later; Japp = — 2.58, D = 5 X 
10~^). The scale in the time series is 250. A solid/dashed 
red line represents a stable/unstable fixed point. Solid/open 
circles represent stable/unstable limit cycles. There is another 
stable steady state at ~ — 8. 



large amplitude relaxation oscillations. Depending on the 
value of the control parameter b, the deterministic sys- 
tem {D = 0) takes one of three solutions: a stable steady 
state for 5 < 5h ~ 0.31535 where 6h is the Hopf bifur- 
cation value, stable small amplitude oscillations (STOs) 
for 6h < < 6c « 0.31854 with be the canard value 
(see, e.g., Ref. HI), and large amplitude oscillations for 
values above the canard value b > b^. As discussed in 
Ref. 26, applying noise to this system {D ^ 0) can pro- 
duce MMOs as noise drives the system between small and 
large amplitude oscillations. 

The infiuence of noise on the MMO dynamics of the 
vdP model has been studied recently in detail^i, pro- 
viding scaling relationships between the proximity to the 
Hopf value and the noise levels that lead to different types 
of behavior. Here we use the model given by Eqs.|4]and[5] 
as a basis for comparison with the biophysical model for 
two reasons. First, the use of g{x) gives control over 
the time scale and shape of the spike (see above). The 
other important effect that the nonlinearity has on the 
system is a shift of the canard point, enlarging the pa- 
rameter region of stable STOs (i.e., bc — bn) by a factor of 
roughly 1.5. This leads to MMOs with well-defined STOs 
over significant parameter ranges even in the presence of 
stronger noise values, in contrast to the vdP model stud- 
ied in Ref. 27. One other difference, compared to the 
vdP model, is that the choice oi g{x) yields a less abrupt 
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transition from STOs to relaxation oscillations, allowing 
for a greater probability of observing spikes with reduced 
amplitude, an effect rarely observed in the SC model. 

As in Ref. [H we augment the MFN model (Eqs. H 
and [S]) with a slow variation in the control parameter 
given by an additional dynamical equation for b. The 
presence of a slowly varying control parameter is a feature 
common to a number of MMO models including the SC 
model, and is necessary to reproduce some of the features 
observed in IF models. Additional comments about the 
model of Eqs. |4] and [5] relative to the augmented models 
are included in Subsec. IV Al 

The first augmented model we consider is an IF model 
with a linear ramp for b plus reset at a fixed threshold 
value for u: 



b = €2- 



(6) 



Once the trajectory surpasses uth, b is reset to 6rs and 
we choose to reset u and v such that the system is on the 
nuUcline (urs = &rs and v^s = &rs(&rs — — &rs)) or very 
close to it (see Subsec. IV Al for a discussion of the dynam- 
ics when reset values are off of the nuUcline). This sys- 
tem, which we reference as the linearly augmented mod- 
ified FN model (LMFN), generates only MMOs within 
the parameter values considered. The stochastic LMFN 
model exhibits features of both a SPHB and CR, depend- 
ing on the value of 62 as well as where the trajectory is 
reset. This allows for a comparison with simple IF mod- 
els, highlighting the effect of noise on the dynamics for 
different ramp speeds of the control parameter b and reset 
values. 

The second augmentation of the MFN model consid- 
ered is a nonlinear it-dependent equation for b{u) with a 
functional form similar to that in the SC model: 



b = 



1 



cxp 



( u-Q.2 \ 
\ 0.1 / 



Cbb 



5 exp 



0.8 



0.15 



(7) 

This augmentation allows for a comparison of cases where 
the dependence of the control variable on u results in 
qualitatively different stable behaviors for the determin- 
istic system. Fig. [2] shows the three possible long time 
deterministic behaviors (D — 0) depending on the pa- 
rameter Cf), keeping the other constants in Eq. [7] fixed. 
For Cb > Cb^H ~ 1-53 there is a stable steady state, for 
intermediate values 1.04 < Cf, < 1.53 there are sustained 
STOs, and for Cb < 1.04 the attracting state is MMOs. 
The number of STO periods between spikes depends on 
Cb and on the other constants in Eqs. HI [5l and [T] We 
refer to this augmentation of the MFN model as nonlin- 
early augmented modified FN (NLMFN) and in the noisy 
version of this model similar patterns as in the SC model 
can be observed. For Cb > Cb^ti noise can excite CR-type 
STOs and for Cb ^ 1.04 the trajectory passes through the 
underlying HB and canard transition thereby displaying 
features of the SPHB. 
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FIG. 2. The bifurcation diagram of the MFN model with tra- 
jectories of the deterministic NLMFN model superimposed. 
The three figures show the three different dynamics of the 
deterministic model depending on the value of cj, (see text): 
top: ct — 1.54, stable steady state; middle: ci, — 1.5, stable 
STOs; bottom: Cb = 1.03, stable MMOs. The timescale in 
the time series is 2. Lines and symbols as of Fig. [1] 



C. Measures 

For the analysis of the MMOs and specifically the 
STOs in the ISIs, we use a number of measures: ISI 
density, amplitude of the STOs preceding a spike, power 
spectral distribution (PSD) of the STOs, and a coherence 
measure. The time series were generated by numerically 
integrating the coupled differential equations given in the 
preceeding subsections using a simple Euler-forward rou- 
tine with constant time step. For the output, it was 
ensured that the data was sampled at a high enough 
rate 1/Ai, such that the STOs could clearly be identi- 
fied (40At <TsTo)- 

The density of the interspike intervals (ISI) was ob- 
tained in the following way: In IF-type models (3DSC 
and LMFN), ISI duration is measured from the reset to 
the crossing of the threshold corresponding to initiation 
of the spike. The typically long refractory period with- 
out STOs in the 3DSC model was removed. In spiking 
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models (7DSC, NLMFN), it is measured as the time be- 
tween two consecutive crossings of a threshold (i.e., the 
time between two spikes including a spike). 

To compare the different models (with different time 
scales), we measure the ISI duration in Tsto, the ap- 
proximate mean STO period. For simplicity, we use only 
one TgTO for each of the model versions and obtain these 
approximately from one time series. We use the follow- 
ing approximate values for Tstq: 3DSC: 107.5, NLMFN: 
0.47, LMFN: 0.45, but note that the actual average Tsto 
vary by about ±5% among the classes introduced in the 
next section. Average ISI duration is used as a coarse 
measurement to calibrate FN-type models to the bio- 
physical (SC) model, and the ISI density is used to com- 
pare models and parameter choices, including influence 
of the noise. 

A second characteristic of the STOs used to calibrate 
FN-type models to the biophysical model is average STO 
amplitude trend before a spike. For this, the maxima 
of the low-pass filtered time series were extracted and 
averaged over many ISIs (for details see App. [C|) . 

To obtain the power spectral distribution (PSD) of the 
STOs, the spikes were removed from the time series and 
it was ensured that the normalization of the PSD is con- 
sistent within each model. From the PSD, we compute a 
measure of coherence, namelyii 

/3 = /^P^ (8) 

where fp is the frequency of the peak in the PSD cor- 
responding to the STOs, hp is its associated maximal 
power and A/ is its full width at half maximum. Since 
the different models use different units, a comparison of 
the absolute value of /3 across models was not performed. 
For details regarding the computation of the PSD and /3, 
see App. [C] 



III. ANALYSIS OF MMO CHARACTERISTICS: 
INTER-MODEL COMPARISON 

In this section we present an analysis of the charac- 
teristics of MMOs with computational measures, focus- 
ing on the STOs, produced by the different models as 
introduced in Section |lll STO amplitude trend before 
a spike and average ISI provide a coarse classification of 
the MMOs in the subsections below. For each class of be- 
havior the parameters of the 3DSC model are chosen to 
produce a type of STO observed in the 7DSC model, and 
the parameters for the FN-type models are chosen such 
that they reproduce the mean ISI and amplitude trend 
similar to the 3DSC. The results illustrate features of the 
MMOs that can be reproduced by both the biophysical 
and phenomenological models. The measures introduced 
in Subsec. Ill CI provide a finer comparison of the three 
broad classes. We introduce different noise levels to high- 
light the mathematical elements of both the underlying 
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FIG. 3. Class 1 MMOs in the LMFN model: trajectory in 
the underlying 2D bifurcation diagram and time series. £2 ~ 
0.0147, D = 10"*, reset to b^s ^ u ^ 0.315, u = -0.12603. 
The scale in the time series is 2. 



deterministic model and stochasticity that influence the 
MMO characteristics. 



A. Class 1: Short ISI, Increasing STO amplitude 

In Class 1, we analyze MMOs with short ISIs during 
which the amplitude of the STOs increases. Examples 
for this class of MMOs are shown in Fig. [T] (top panel) 
for the 3DSC model and Fig. [3] for the LMFN model. 
The comparison of the amplitude dynamics of the STO 
and ISI density is shown in Fig. H for 3DSC and LMFN, 
where the ISI density and mean agree up to a shift due 
to different reset conditions. This Class 1 behavior is not 
readily reproducible within the NLMFN model where the 
choices of Cb needed for short ISIs lead to STOs with a 
larger amplitude throughout the ISI. 

As shown in Sec. |IT1 the underlying dynamics for both 
the 3DSC and the LMFN model corresponds to a slow 
passage through a Hopf bifurcation. In the SC model, the 
increasing amplitude of the STOs are generated when the 
system spirals away from an unstable fixed point for Vs 
values beyond the subcritical HB. In the LMFN model, 
the trajectory passes through the underlying supercriti- 
cal HB at &H followed by an increasing STO amplitude, 
with escape from STOs to a spike appearing for relatively 
large b (around 0.34) beyond the canard point be- The 
significance of the speed of variation of b in the LMFN 
model is discussed further in Subsec. IV Al 

The noise levels shown in Figs. [1] (top panel) and [3] 
are among the lowest at which stochastic effects are ob- 
served, with the STO dynamics dominated by the de- 
terministic trend of increasing amplitude. Increasing the 
noise strengths has very similar effects in both models. 

• Typical for a SPHB, increased noise drives an earlier 
escape, indicated by shorter ISIs (Fig. IH) and lower av- 
erage amplitude before escape (Fig. [5|) . Differences due 
to the nature of the underlying Hopf bifurcation (sub- 
/supercritical) are evident only in the amount of reduc- 
tion of the amplitude of the STOs. This reduction is 
more significant in the 3DSC, where the STOs are un- 
stable for the underlying 2D system with > r^.n, as 
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FIG. 4. Top: Average amplitude before a spike (obtained 
by finding the maxima (see Subsec. Ill C|) ) for time series of 
Class 1. A (solid red line): 3DSC (left ordinate); o (dotted 
black line): LMFN (right ordinate). Bottom: Densities of 
ISI lengths for long (~ 5500 spikes) versions of the time series 
shown in Figs. [1] (top panel) and|3] Solid line (red): 3DSC; 
dashed line (blue): LMFN. Parameters: lapp = -2.45, D = 



10"^ (3DSC); £2 = 0.0147, D = 
0.315, w = -0.12603 (LMFN). 



10" 



reset to &rs — 



compared to the LMFN model, where the supercritical 
Hopf bifurcation yields stable STOs for b > bn. 

• Increased noise shifts the ISI density to include peaks 
at shorter durations, consistent with earlier escape de- 
scribed above. Contributions at longer ISIs are due to a 
small fraction of trajectories for which the noise disrupts 
the transition and lengthens the ISI. Class 1 behavior 
survives only for low to moderate noise levels, as larger 
noise regularly eliminates the STOs. This is illustrated 
by the considerable spread of both the ISI density and 
the PSD (Fig. [7]) for large noise levels. 

• The coherence measure /? is not clearly defined for 
MMOs of Class 1, since the strong amplitude trend of 
the deterministic dynamics within the short ISI periods 
yields PSDs with many well-defined peaks (one of them 
at the STO frequency - see Fig. [7]). The SC model shows 
more peaks due to larger amplitudes of the STOs across 
the ISI. The strong peak at low frequencies comes from 
the spiking (i.e., thresholding and concatenation). 



B. Class2: Long ISI, increasing STO amplitude 

Next we analyze STOs appearing in longer ISIs with 
a characteristic increasing amplitude before a spike. We 
consider two examples from the 3DSC model that have 
similar amplitude trends but different ISI behavior in the 
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FIG. 5. Average amplitude before a spike for various values 
of D. Top: 3DSC; bottom: LMFN. The parameters are as 
of Fig. H 
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FIG. 6. Densities of ISI lengths for various values of the noise 
strength D. Top: 3DSC; bottom: LMFN. The parameters 
are as of Fig. [H 



presence of noise. As discussed in Sec. |TT] the choice of 
^app for these two examples, /app = —2.55 and /app — 
—2.58, corresponds to different underlying deterministic 
behavior, regular MMOs and stable steady state, respec- 
tively. We refer to these as Class 2 A and Class 2B, re- 
spectively. 

Sections of the time series as well as the trajectories 
within the underlying bifurcation are shown in Figs. [T] 
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FIG. 7. Power spectral density (PSD) of the time series in 
Class 1 for various values of D. Top: 3DSC (/app = —2.45); 
bottom: LMFN (e2 = 0.0147, reset to brs = u = 0.315, v = 
—0.12603). For details of how we computed and normalized 
the PSD see App. O 



FIG. 8. Examples of the dynamics of the MFN model for 
Class 2. Top: NLMFN with Cb = 1.1, D = 5x 10"*; bottom: 
LMFN with £2 = 0.00055, D = 3 x 10"^ reset to 6rs = 0.315. 
Trajectory in the underlying 2D bifurcation diagram and time 
series. The scale in the time series is 5. 



(bottom panel) and M for the 3DSC, NLMFN, and 
LMFN, respectively. These time series were calibrated 
for similar average amplitude trend and average ISI du- 
rations. Fig. ini shows that for low noise levels, it is possi- 
ble to find parameter values for LMFN that compare well 
in amplitude behavior with 3DSC for both Class 2A and 
Class 2B. For NLMFN, the amplitude behavior is similar 
to 3DSC for Class 2B but differs in Class 2A for STOs 
well before the spike. In the refractory period following 
the spike in the NLMFN the system relaxes near but not 
on the left null cline with 5 < 6h , yielding slowly damped 
STOs in the first part of the ISI. For the LMFN, the re- 
set value for u and v is chosen very close to the fixed 
point, avoiding this part of STOs with decaying ampli- 
tude. Comparable average ISIs are easily obtained in the 
LMFN with €2 chosen an order of magnitude smaller than 
in Class 1, and for NLMFN with Cf, near unity for weak 
noise. Differences in the ISI density are apparent from 
Fig. [TU] and are related to differences in the deterministic 
dynamics as discussed in Sees. llVl andlVl 

As we see below, by tuning parameters and noise lev- 
els, the FN-type model can capture some but not all of 
the amplitude, ISI and coherence behavior for Class 2. 
Results for different noise levels point to differences in 
the underlying structure of the models. 

• For noise levels increased by an order of magnitude, 
both subclasses for 3DSC show roughly a 20-30% de- 
crease in the average amplitude of the largest STO imme- 
diately preceding a spike, while the FN-type models show 
a slightly greater decrease (Fig. [TT]) . For Class 2B, the 
entire average amplitude curve of 3DSC shifts to lower 
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FIG. 9. Average amplitude before a spike for time series from 
Class 2 from three different models: A (solid red line): 3DSC 
(~ 1400 spikes); x (dashed blue line): NLMFN (~ 2500 
spikes); o (dotted black hue): LMFN (~ 1000 spikes). Top: 
Class 2A: 3DSC: /app = -2.55, D = 10"^; NLMFN: ct = 
0.95, D = 10"**; LMFN: ea = 0.001, D = 10-^ 6rs = 0.3155; 
bottom: Class 2B: /app = -2.58, D = 5 x 10"^ (Fig. [1] 
bottom panel); NLMFN: a = 1.1, /) = 5 x 10"* (Fig. [8] 
top panel); LMFN: £2 = 0.00055, D = 3x 10~^ = 0.315 
(Fig. [H bottom panel). 
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FIG. 10. Densities of ISI lengths for Class 2 MMOs. Solid line 
(red): 3DSC; dashed line (blue): NLMFN; dotted line (black): 
LMFN. Top: Class 2A; bottom: Class 2B. For parameters, 
see Fig. (9] 



levels with increased noise, while both FN-type models 
show both shifts and changes in the average amplitude 
curve shape. For Class 2A, both 3DSC and LMFN show 
an increase of the average amplitude of STOs well be- 
fore the spike for increased noise, indicating STOs driven 
both by CR and by SPHB. For the NLMFN model, the 
return sets b < bn following an earlier spike, so that the 
combination of damped STOs for b < bn and CR effects 
yields robust STOs with less variation with noise level 
than for for the LMFN model. 

• In Class 2B, the 3DSC model has longer tails in the 
ISI density for both small and intermediate noise levels. 
These longer tails are not seen in the FN-type models due 
to differences in the underlying deterministic dynamics 
(see Sec. lIVp . For weak noise, the multi-modal ISI den- 
sity for the NLMFN is due to the stronger attraction 
to the underlying deterministic MMO or STO dynam- 
ics, with larger STOs initially in the ISI. In the LMFN 
model the reset onto the nullcline avoids damped STOs 
initially in the ISI, allowing a unimodal ISI density. For 
both Class 2 A and 2B, 3DSC and LMFN show a shift 
to shorter ISIs and a similar shape of the ISI density for 
larger noise, with a larger variance for LMFN in Class 
2A. For NLMFN there is no shift in the ISI density with 
increased noise, just an increased variance both for Class 
2 A and 2B, reflecting the dynamic return mechanism fol- 
lowing the spike. 
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FIG. 11. Average amplitude before a spike in the two subclasses of Class 2 for various values of D. Left column: 3DSC; 
middle column: NLMFN; right column: LMFN. Top row: Class 2A; bottom row: Class 2B. For parameters, see Fig. |9] 
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FIG. 12. Densities of ISI lengths in the two subclasses of Class 2 for various values of noise parameter D. Left column: 3DSC; 
middle column: NLMFN; right column: LMFN. Top row: Class 2A; bottom row: Class 2B. For parameters, see Fig. |§] 
In the top row, middle panel, values for large p are cut off. 



• Fig. [13] shows the differences between the coherence 
measure /3 (Subsec. Ill Cp for the three different models 
for various noise levels. For 3DSC in Class 2B the STOs 
are driven by CR and there is a peak in /3 (indicating 
CR) , while for Class 2 A noise reduces (3 in general, since 
the STOs are driven by a SPHB. For low noise levels in 
NLMFN, the robust STOs are represented by a multi- 
peaked PSD (equivalent to the ISI density - see above). 



for which /? is not defined. For larger noise levels in Class 
2, /3 is defined, but noise only disrupts the underlying 
STOs, so (3 decreases with noise. For LMFN, a mod- 
erate noise level can enhance the STOs, particularly at 
the beginning of an ISI, yielding an optimal noise level 
for coherence in both types of Class 2. The influence 
of the underlying bifurcation structure on the coherence 
measure is discussed further in Sec. IIVI 
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FIG. 13. Coherence measure /3 (see Subsec FlI C[) for the 3DSC 
model (top), NLMFN (middle) and the LMFN model (bot- 
tom) for various values of noise strength D. Class 2A is 
within the red solid lines, Class 2B within the blue dashed 
lines, and Class 3 within the black dotted lines (blue dashed 
for NLMFN). 



FIG. 14. Trajectories of Class 3 MMOs generated with the 
following models: top: 3DSC (/app = -2.7, D = 10""); 
middle: NLMFN (c^ = 1.1, D = 10"^); bottom: LMFN 
(£2 = 0.00055, D = 10~^ 6rs = 0.305). The scales in the time 
series are 250 (SC) and 5 (MFN). 



C. Class 3: Long ISI, constant average STO amplitude 

Here we analyze STOs with an average ISI similar in 
length to Class 2, but with a constant average amplitude 
preceding a spike. Fig. [Til shows the trajectories from the 
3DSC, NLMFN, and LMFN calibrated for similar aver- 
age amplitude behavior and average ISI duration. Class 
3 can be observed in 3DSC only at higher noise levels 
compared with Classes 1 and 2. 

Fig. [15] shows that the average STO amplitude for 
Class 3 is constant for at least the last ten periods of 
STOs with a slight amplitude increase in transition to 
the spike. This nearly constant average amplitude re- 
sults from strong variation in amplitude with a constant 
ensemble average. A clear single peak in the PSD (data 
not shown) confirms well-defined STOs with variability 
in phase and amplitude. 

The amplitude and average ISI characteristics of these 
MMOs are readily reproduced with the FN-type models 



at stronger noise levels with the average amplitude lower 
by roughly a factor of two compared to Classes 1 and 2. 
Class 3 behavior can also be obtained for the MFN model 
(Eqs^^ and [5]) with constant control parameter h as in 
Ref. [2y. The ISI densities shown in Fig. [TH] are similar 
to each other, with a slightly more significant tail in the 
ISIs for 3DSC. 

Figs. [TBI and [Tfl show the trend of the STO amplitude 
and the density of ISI when noise levels increase. 

• For stronger noise the increase in amplitude directly 
before the spike in the FN-type models is eliminated, 
with stronger noise levels rather than deterministic dy- 
namics dominating the transition to spikes. At higher 
noise levels, distinguishing stochastic fluctuations from 
regular STOs is less dependable, so that the amplitude 
dynamics for the largest noise levels are not shown in 
Fig. [m In the 3DSC model, there is a slight increase 
in average amplitude during the ISI which is due to our 
algorithm (Subsec. Ill C[) not fully eliminating the trend 
of increasing V in steady state with increasing for 
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FIG. 15. Top: Average amplitude before a spike for the time 
series from Class 3: A (red solid line): 3DSC (~ 1600 spikes); 
X (blue dashed line): NLMFN (~ 1200 spikes); o (black dot- 
ted line): LMFN (~ 1200 spikes). Bottom: Densities of ISI 
lengths for Class 3 MMOs. Solid line (red): 3DSC; dashed line 
(blue): NLMFN; dotted line (black): LMFN. For parameters, 
see Fig. [TH 



• The ISI densities for 3DSC and LMFN are more con- 
centrated at shorter ISI durations, while the ISI density 
for NLMFN remains spread over a large range of ISIs. 
This is due to the underlying stable STOs for NLMFN, 
together with the nonlinear return mechanism that typi- 
cally returns b to values well below 6h following a spike, 
thus allowing for a range of ISI durations. With the 
ISI density concentrated at low values, clustered spikes 
are more frequent in 3DSC and LMFN with nearly tonic 
spiking at the higher noise level (cf. Ref . . 

• For Class 3, /3 for the 3DSC model shows a similar 
behavior as for Class 2B (Fig. [13]), with a clear maximum. 
The values of /3 in Class 3 are lower than in Class 2 for 
3DSC, mainly because of smaller amplitudes of the STOs 
(being inversely proportional to the distance from rs^u for 
CR—- - also see Subsec. IIIC|) . Similarly, the trend for /3 
in Class 2B and 3 is the same in LMFN. For the NLMFN 
model, the value of Cb is the same for Class 2B and Class 
3 and the behavior of /3 is shown in Fig. [131 



IV. INTRA-MODEL COMPARISONS FOR STELLATE 
CELLS 

The comparisons in the previous section point to the 
role of the underlying bifurcation structure on the charac- 
teristics of the STOs and MMOs. We summarize the re- 
sults for the 3DSC model from Sec lIIIi comparing within 
that model. We also discuss alternative values of /app for 
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FIG. 16. Average amplitude before a spike in Class 3 for 
various values of D. Top: 3DSC (/app ~ —2.7); middle: 
NLMFN {cb = 1.1); bottom: LMFN (e2 = 0.00055, 6rs = 
0.305). 



the 7DSC model in order to compare with previous stud- 
ies of the full model. 



A. Comparison for different ranges of /app 

Within the 3DSC model (and similarly the 7DSC 
model) differences in the stochastic behavior can be re- 
lated to the differences in the deterministic behavior for 
different ranges of Japp. For /app < /app,H « -2.575 



in 3DSC (see Subsec. fll Ap . the stable steady state has a 
value of = rs^o well below Tg^. For /app > /app.H there 
is no steady state for the deterministic system, and the 
STOs are driven by SPHB with as the slowly varying 

char act er- 



app 



^app,H 



control parameter. In the case /; 
istics of both CR and SPHB are observed. 

• Amplitude: For /app well below /app.H, Ts,o is well 
below Ts^K, so only high noise levels can drive MMOs, 
and STOs are purely of the CR-type without any trend 
in the average amplitude (as in Class 3). In contrast in 
Class 1, with /app clearly above /app,H the STOs due to 
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FIG. 18. ISI densities from time series generated with the 
7D version of the SC model (using the alternative equation 
(Eq. IB9[) 1 with various values of Japp (—2.435 in the top 
panel corresponding to family Jo(Ji)^; —2.4416 in the bot- 
tom panel corresponding to family Jo(Ji)*). The determin- 
istic case {D — 0) is shown as a histogram (solid red lines) 
whereas the data for simulations with noise {D ^ 0) is shown 
as line plots (dashed blue/dotted black lines). We adopt the 
notation for the MMOs from Ref. Tsto « 99. 



FIG. 17. Densities of ISI lengths in Class 3 for various values 
of the noise strength D. Top: 3DSC; middle: NLMFN; 
bottom: LMFN. For parameters, see Fig. 1161 



the SPHB have an increasing trend in amplitude before 
the spike. Well-defined periods of STOs survive only for 
lower noise levels, as the slow passage is sensitive to very 
weak noise. For values of /app closer to /app,H, as in 
Class 2, both CR and SPHB can affect the dynamics, 
so that periods of STOs with increasing amplitude trend 
can survive over a large range of noise levels. 

• ISI density: Significant differences are seen at lower 
noise levels, where longer tails in the ISI density corre- 
spond to STOs of the CR-type, not observed in STOs 
driven by SPHB. 

• Coherence measure /3: CR-driven STOs in Class 2B 
and Class 3 are responsible for the more pronounced peak 
in the coherence measure /? in Fig. 1131 

Intermediate values of /app = —2.5, —2.65 were also 
analyzed (not shown), and as expected intermediate be- 
havior between the three classes is observed. For /app — 
—2.5, STO and ISI behavior was closer to that of Class 2B 
than Class 1, without contributions to longer ISIs, and 
for intermediate to larger noise levels very weak coher- 
ence is observed with a mild peak in /? around D sa 10~^. 



B. Families of MIVlOs for Japp > /app.n at weak noise 

Here we consider the effect of noise on MMO families 
in the full 7DSC model in the range /app > /J°p jj (i.e., 
within the regime of deterministic MMOs), as illustrated 
through the ISl density. The 3DSC model shows a similar 
behavior, with a shift in the values of /app- Recall that 
for any fixed /app > /J^p^H ~ —2.702, the underlying de- 
terministic dynamics is a slow passage of Ts through r-g^. 
We restrict our study to very low noise levels and con- 
sider trajectories that fit into Classes 1 and 2A (Sec. IIII[) 
when stronger noise is added. For the sake of compar- 
ison with the analysis in Ref. i34l . we use an alternative 
equation for the dynamics of r^. Instead of Eq. [31 we 
use Eq. IB9I given in App. [B] which alters the numerical 
values of the model but not the qualitative features. 

Fig- EH shows examples of the ISI density for different 
families of MMOs. The 1" families (for notation, see 
Ref. ^) have MMOs with s STOs in each ISI, yielding 
ISl densities with a single point mass. The MMOs for a 
second, more complex family, Ji(Jj)', have an ISI with 
i STOs followed by / ISls each with a different type of 
j STOs that increase in amplitude and period for each 
of the subsequent I ISIs. The Ji[,Jjf families then have 
1 + / point masses in the ISl density. We consider only 
those MMO families that are stable in the deterministic 
model. In the range —2.55 < /app < —2.4 stable 1^ 
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FIG. 19. As of Fig. [18] with J^pp = -2.5 (deterministic fam- 
ily 1^ - top) and /app — —2.66 (deterministic family 1^^ - 
bottom). 



and JiiJjY families are found for alternating ranges of 
lapp- For /app < —2.55, the stable solutions are only 
1" families densely filling the space in /app'^^. The time 
series analyzed earlier as Class 1 in this notation 'would 
be a mixture of 1^ and 1'^. 

For values of /app near the center of the stability 
range for a particular type of MMO, 'well separated from 
other families, -weaker noise broadens the ISl density 
and stronger noise shifts the density to-wards shorter ISls 
"when noise drives an early escape typical in a slo"w pas- 
sage through a HB, as sho'wn in Fig. [18] for Jo(Ji)^. The 
bottom panel of Fig.[T5]sho'ws that for values of /app near 
the edge of a stability range, very "weak noise can lead to 
a mixture of close families of MMOs. As noise drives a 
faster escape to spiking the peak corresponding to the 
longest ISI and largest STO of Jo(Ji)'* is reduced as the 
nearby MMOs of Jo{Ji)^ occur, appearing as stronger 
nearby peaks in the ISl density. There is also consider- 
able sensitivity in the ISI density "when stability ranges 
are small and densely packed. Noise can drive a sampling 
of a mixture of families close in terms of the bifurcation 
parameter, yielding clearly defined peaks that are not 
present in the deterministic case, as sho'wn in Fig. [TO] In 
the top panel, the deterministic MMO is 1^, and noise 
excites a second peak for a family 1^ or a mixture of 
families bet'ween 1^ and 1^. At /app = —2.66 the stable 
solution is 1^^ and 'weak noise drives a sampling of the six 
nearby MMO solutions 1* 'with 10<s<15. At stronger 
noise the longer ISls are lost and the ISl density spreads 
out and shifts to'wards shorter ISls. 

In conclusion, even very "weak noise can alter the dy- 
namics of a MMO-generating system significantly, partic- 



ularly if the deterministic solutions are close in parame- 
ter space. Weak noise drives a sampling bet"ween a fe"w 
or many of these solutions. Then it can be difficult to 
distinguish bet'ween complex deterministic solutions like 
the Ji{JjY families and a noisy trajectory sampling a fe'w 
different MMO families. 



V. INTRA-MODEL COMPARISONS FOR THE 
MODIFIED FN MODEL 

As in the SC model, the underlying bifurcation struc- 
ture and deterministic behavior influence the source and 
characteristics of the STOs. For the FN-type models, the 
underlying bifurcation is a supercritical Hopf, ■with sta- 
ble STOs in the 2D reduced system for u-v 'with constant 
6h < 6 < be- 



A. Linearly augmented modified FN model 

In the LMFN model 'with appropriate reset the slo'wly 
varying control variable S"weeps through an underlying 
HB and canard transition, "which al'ways yields MMOs 
in the deterministic setting. We consider t"wo aspects 
of the control parameter that can be varied to change 
the characteristics of the STOs, the rate of change in 
b, b — €2, and its reset b^s, foUo'wing a spike/threshold 
crossing. 

In the different classes studied in Sec. Illli we analyzed 
four different combinations of 62 and bj^s'- slo"wer variation 
in b given by smaller £2 — 0.00055, "with reset 6rs < 
(Class 3) or b^s ~ 6h (Class 2B), and larger values of 
£2 = 0.001 (Class 2A) and £3 = 0.0147 (Class 1) both 
with b^s ~ bn. 



1. Characteristics for £2, 6rs in Section [TTTI 

The main underlying feature that influences amplitude 
dynamics for these different parameter combinations, is 
time spent in the parameter range b > b^. • Ampli- 
tude increase is more pronounced for low noise levels and 
smaller values of £2. For low noise levels, the system 
follows the underlying STO dynamics, less likely to be 
driven to spiking by noise. For smaller values of £2, the 
control parameter behaves almost as a constant relative 
to the oscillation frequency. This allows attraction to 
larger STOs for slowly increasing 6h < b < be, with the 
possibility of reaching values of 6 > 6c before spiking, 
so that larger amplitudes are typical before escape. For 
larger values of £2, as in Class 1 or Class 2 A, b increases 
through the STO region less slowly. Then the trend of 
increasing amplitude consistent for 6h < ^ < ^'c survives 
but increased variation in b limits the attraction to the 
stable STOs as b approaches and exceed be- The behavior 
is closer to a series of transients that is more susceptible 
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to noise-induced escape to spiking. For larger noise lev- 
els, the escape to spiking is noise-driven in or near the 
stable STO region. An increase in STO amplitude far- 
ther from the spike, for 6 « 6h is observed due to the role 
of CR together with the stable STOs. An earlier escape 
limits the opportunity for amplitude increase closer to 
the spike. 

• The relation between £2 and &rs and ISI density is 
as expected: smaller values of £2 or &rs lead to greater 
likelihood of longer ISIs. 

• The coherence measure (3 shown in Fig. [13] increases 
for long stretches of large amplitude oscillations, and de- 
creases with fluctuating amplitude. Trajectories with the 
reset close to the HB 61s ~ &h show the largest value of fi 
for moderate noise levels. Slower variation in b in these 
cases also reduces fluctuation in amplitude, thus increas- 
ing coherence. For reset h^s < &h there is fluctuation 
in amplitude due to oscillations driven by both CR and 
STOs, so /3 is reduced somewhat overall. For larger val- 
ues of noise coherence is destroyed in all cases. 



2. Characteristics for slower increase in b with bis < bn 

Fig. [20] shows the analysis of an additional parame- 
ter combination, a very slow sweep through the HB and 
reset before the HB (e2 = 0.0001, 6rs = 0.305), whose 
coherence measure /3 was included in Fig. 1131 

• The difference in amplitude dynamics for lower and 
higher noise levels is as described above. For lower noise 
levels and very slow dynamics of b, the system is attracted 
to larger amplitude STOs before spiking, so that larger 
amplitudes are typical for more STOs before escape. 

• The ISIs are very long, and only for stronger noise are 
they comparable to Class 3 (cf. Fig. [T5)) . In contrast to 
Class 2, only for stronger noise levels is the ISI density 
concentrated at shorter ISIs, since the trajectory takes 
longer to reach &h- 

• The shape of the curve /3 vs. Z? is qualitatively sim- 
ilar to the classes studied in Sec. IIIII However, at very 
low noise values, /? is larger than other cases, caused by 
CR-drivcn STOs and slower variation in b as described 
above. 

The last panel in Fig. [50] shows the ISIs for an interme- 
diate case, with larger £2 and 6rs < &h- The ISI density 
then shows characteristics between Class 2 and Class 3, 
with a narrow ISI density for smaller noise levels, and 
longer tails for stronger noise, as oscillations of the CR- 
type appear for 6 < 6h- 

We add one additional note about resetting off of the 
nuUcline. The choice of reset used in Sec. IIIII (reset very 
close to the nullcline), can be relaxed to a reset near 
the nullcline without qualitative changes in the measures 
used above. A choice of reset away from the nullcline 
modifies the amplitude dynamics of the STO at the be- 
ginning of the ISI. For reset with b^s well below bn, STOs 
with considerable amplitude follow the reset and relax as 
b approaches 6h, similar to NLMFN. For Classes 1 and 
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FIG. 20. Average amplitude (top) and ISI density before a 
spike (middle) for the LMFN model for a very slow sweep 
through the HB (e2 = 0.0001) and reset before the HB (6rs = 
0.305). ISI density (bottom) for an intermediate case (e2 ~ 
0.0055, b,s = 0.305). 



2 this increased amplitude supports an earlier escape to 
spiking and reduced ISIs, particularly in Class 1. For low 
noise levels in Class 2 the ISIs are then narrower, since 
the larger STOs are less susceptible to small noise. For 
Class 3 the stronger noise dominates, so that variation in 
the reset has limited effect. 



As noted in Subsec. lilB| the model of Eqs. [4] and [5] has 
a larger distance between bn and b^ as compared with 
vdP, yielding a larger region in parameter space where 
well-defined periods of STOs are observed in the pres- 
ence of noise. This behavior is confirmed by comparisons 
of PSDs for MFN and for a rescaled version of vdP to 
effectively increase the canard value. This comparison 
illustrates that the distance between 5h and be plays a 
critical role in the robustness of STOs and MMOs, in 
addition to the reset and ramp speed. 
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B. Nonlinearly augmented modified FN model 

As described in Subsec. IIIB[ the long time behav- 
ior of the NLMFN deterministic model includes differ- 
ent states depending on C{,: C}, > 1.53: steady state; 
1-04 < Cb < 1.53: stable STOs; Cb < 1.04: MMOs. 
In Sec. mil we have considered only those values of Cb 
that support phenomenological behavior observed in the 
SC model. Here we also discuss cases that differ qualita- 
tively from those studied in Sec. lIIIi and compare the am- 
plitude, ISI density, and coherence measure within this 
model. 



1. Values of Cfc « 1 

Values of Cfc = 0.95 and Cf, 1.1 were used in Sec. lIIII to 
approximate Class 2 behavior. For Cb = 0.95 the under- 
lying deterministic dynamics is MMOs, with b regularly 
crossing be- Then the observed measures are consistent 
with a slow passage through a Hopf bifurcation and ca- 
nard point as observed in the other models: 

• Increasing amplitude before the spike over a range of 
noise levels. 

• Multi-peaked ISI density for lower noise levels and 
concentrated ISI density for larger noise. 

• Decreasing trend in /3, defined only for larger values 
of the noise. 

For values of Cb < 0.95 the return value of b following a 
spike takes values between 6h and be, thus shortening the 
ISI significantly. As discussed above, we do not classify 
this behavior as Class 1, since the return value yields 
larger amplitude STOs and limited increase for low noise 
levels, similar to those observed for a related model in 
Subsec. rvTBl 

In contrast, for Cb — 1.1, the deterministic system 
exhibits stable STOs where b oscillates near be- The 
stochastic behavior is similar to that of the case with 
Cb = 0.95, but differences due to the underlying deter- 
ministic dynamics are observed for lower noise levels and 
Cb = l.l: 

• Increased noise drives an earlier escape to spiking. 
This results in a reset oi b < 6h, yielding both CR-type 
and SPHB-driven oscillations, and a limited opportunity 
for amplitude increase just before the spike. 

• The ISI density is broader for Cb = 1.1 (compared to 
Cb = 0.95). 

2. Intermediate values: ct = 1.3 

The attracting deterministic behavior for b consists of 
small oscillations midway between 6h and be, so charac- 
teristics include some elements from Class 2 for smaller 
noise, and similarities to Class 3 and Cb = 1.5 (see below) 
for larger noise: 

• The amplitude dynamics is similar to Class 2 for 
lower noise, with a weaker increase. 
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FIG. 21. Average amplitude before a spike (top) and ISI 
density (bottom) for the NLMFN model with cj, — 1.6 at 
which the deterministic system is quiescent. 



• The ISI behavior has characteristics closer to Class 

2. 

• The coherence measure /3 vs. D behaves similarly to 
Cb = 1.5, with larger values of (3 for small noise. 



3. Large values of ct = 1.5, 1.6 

For larger values of Cb, the attracting state of the de- 
terministic system takes values of b near 6h- For Cf, = 1.6 
the long term deterministic behavior is quiescent with 
the steady state corresponding to b just below 6h- For 
Cb — 1.5 the deterministic system has stable STOs, with 
small oscillations in b just above 6h. To drive MMOs 
in either of these cases, noise must be strong enough to 
cause excursions beyond the canard transition at be, cer- 
tainly stronger than in the cases for Cf, « 1 above. ISI 
and amplitude measures can therefore only be derived 
for these stronger values of noise, leading only to Class 3 
MMOs with the following characteristics (Fig. [2T|) : 

• Primarily constant trend in average amplitude before 
a spike. 

• Broad ISI density with long tails. 

• A clear difference between Cb — 1.5 and 1.6 due to 
the underlying deterministic dynamics is seen in the co- 
herence measure /3 (Fig.[T31). For Cb — 1.5 all noise levels 
disrupt the coherence of the STOs produced deterministi- 
cally, so that /3 strictly decreases with noise. For Cb = 1.6, 
/3 has a maximum indicating an optimal noise level typ- 
ical of CR-driven oscillations^!^^. 



17 



VI. APPLICATION TO OTHER MODELS 
A. The 3DSC model with noise in /app 

The stochastic nature of the SC system results from 
noise in ion channels and noisy external signals. The 
model we consider above follows Ref. IH, including the 
primary noise source in the gating variable of the persis- 
tent sodium channel. Here, we compare our results with 
those obtained with a different noise source in the 3DSC 
model, namely, in the applied current /app- The con- 
sideration of this additive noise is motivated by several 
factors. First, systems with multiple time scales often 
can be divided into regions of time or space where noise 
plays a more or less significant role. This raises ques- 
tions about where or when noise plays a critical role, and 
whether the type of noise can make a significant differ- 
ence. Second, since /app is the only variable in the system 
that can be easily experimentally controlled, this raises 
the question of whether varying /app can be used as a con- 
trollable input to probe the dynamics. Providing noisy 
input is a technique used in experimental neuroscience 
(e.g., Ref. HI) to investigate the structure of the com- 
plex dynamics. Finally, understanding whether different 
types of noise sources have similar or different effects con- 
tributes both to the understanding of the mechanism for 
the dynamics and to model identification. 

To test whether fluctuations in /app can reproduce dy- 
namics similar to those generated by noise in the persis- 
tent sodium current, we consider the voltage equation in 
the interval near the end of the ISI, when the dynamics of 
V is slow. This is the stage at which STOs are generated 
for parameters near the Hopf bifurcation of the reduced 
system. In that part of the cycle, V can be approxi- 
mated by a constant Vq in Eq. [TJ This approximation 
suggests that even though the noise is parametric, at the 
end of the ISI when the STOs are generated, the noise 
behaves as additive noise to leading order with a coef- 
ficient D' w (0.15Gp(Vb - E^^)fD and Eq. [J can be 
approximated by 



app 



Gl{V^-El)-Gp 



1 



1 



exp(-™ 



X {Va - ^Na) - G,,(0.65ry -f 0.35r,)(l/o - Eh) 
1 



C 



(9) 



Eq.|9]can be interpreted as having a noisy applied current 
with /app -I- \/2D'r]{t) and D' w 68L>'(for Vb w -55) 
instead of the original noisy gating variable. Then one 
would expect that the effects of different noise levels on 
the STOs as described in Sec. HIT] could be reproduced 
by varying the noise level in /app, as long as the results 
are scaled in a manner consistent with Eq. [9] Indeed 
for the ranges of noise levels considered in this paper, 
the results for stochastic /app are essentially the same 
as described in the classification of Sec. IIIII With the 



above scaling, the plots for the coherence measure /? (as 
in the top panel of Fig. [T3]) overlap almost perfectly. We 
also reproduced similar behaviors of the ISI density and 
amplitude dynamics as analyzed for the 3DSC model in 
Sec. HIT] (data not shown). 



B. A self-coupled FN-type model 

Here, we analyze another MMO-generating model, a 
self-coupled modified FN model derived and used for the 
study of coupled Hodgkin-Huxley neuron o^'^'^ . Similar to 
the MEN model presented earlier, v represents a voltage 
variable, h is gating variable, and s is a dynamic coupling 
between neurons. Compared to Refs. [tI andlssl we add a 
noise term in the first differential equation: 



w = - 0.5(t;^ - V + l) + h- jsv + 

h =eDi2h + 2.6v) 

s =l3H{v){l - s) - eoSs. 



(10) 

(11) 
(12) 



H{v) is the Heaviside function. This system of differ- 
ential equations has multiple time scales similar to the 
models analyzed so far: the dynamic variable s plays the 
role of a slowly varying control parameter, decreasing 
and passing through a HB at sh for the underlying v-h 
system. Also, the spikes are included in the model (not 
just a threshold and reset) but this is done using a sharp 
switch {H{v)) rather than a smooth nonlinear function 
for the control variable as in our NLMFN. The param- 
eters used in Ref. are as follows: coupling strength 
7 — 0.5, activation rate /3 — 0.035, slow time scale of 
the relaxation oscillations eo = 0.015, and d = 0.565 
regulating the speed of the dynamic coupling. Ref. an- 
alyzes the deterministic {D = 0) dynamics of the system 
in detail. By varying the parameters in Eqs. llOf - 1121 in 
particular S and D, we find MMOs according to the clas- 
sification scheme of Sec. IIIII We choose en — 0.01 and 
7 — 0.9 to obtain a sharper canard transition, similar to 
the MFN models analyzed earlier. There are a number 
of differences between Eqs. [TUHT^ and the MFN model. 
The slow variation of s is nonlinear, slowing down signif- 
icantly in the SPHB, a behavior somewhat similar to b 
in the NLMFN case with Cb ~ 1.5. Another noticeable 
difference between Eqs. [TUHT^ and the MFN models is 
that the dynamics of s yields a return value well below 
the Hopf point sh following a spike. After the return 
there is a significant interval during which s approaches 
Sh, which has implications for both the amplitude and 
ISI behavior. 

Class 1 time series are obtained with 5 — 2.5 and 
D = 10~^, yielding both the amplitude dynamics and 
the ISI density similar to those shown in Fig. |31 One 
difference for the model of Eqs. [TUl fT^ is that increasing 
noise strength typically increases the amplitude of the 
STOs before a spike. This is consistent both with the 
slowing of s near sh and with the opportunity for CR- 
type STOs in the interval where s > sh- As in Class 
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FIG. 22. The underlying (fixed s) bifurcation diagram of 
the self-coupled FN-type model with e_D = 0.01 and 7 = 0.9 
together with the trajectory and time series for 5 = 0.3 and 
D = 10"^ At lower noise values, the trajectory crosses n = 
at lower values of s (around 0.024 for D = 0). Lines and 
symbols as in Fig. [21 The time scale in the time series is 100. 



1, we see the appearance of new well-defined peaks in 
the ISl density with low noise levels. For this particular 
choice of parameters, the first of these new peaks appears 
at longer ISIs, similar to the behavior in the upper panel 
of Fig. m for the 7DSC model. 

For time series with the characteristics of Class 2, we 
find (5 = 0.3 and D = 3 X 10^^ in Eqs. UnHH yielding 
dynamics comparable to Sec. IIIII in terms of average ISI 
and amplitude trend (Fig. [23l) . For the STOs well be- 
fore the spike, an increase of the amplitude results from 
the algorithm that subtracts the average u as w has a de- 
creasing trend. The amplitude before the spike increases 
with noise level, once again due to the combination of 
CR-driven STOs before s crosses shi and a slowing of s. 
Larger values of D shift the ISI density towards smaller 
values and broadens it, yet a certain minimal ISI is main- 
tained due to the return of s well beyond sh- The effect 
of such a return mechanism was also observed in the ISI 
density in the NLMFN. 

Reduced speed of s (5 = 0.17) combined with stronger 
noise (D = 3 x 10~^) leads to time series that fit into our 
Class 3. The amplitude dynamics is very similar, but the 
ISI density is different from those shown in Fig. (TS] At 
this noise level, the ISI density is still centered around 15 
TsTOi again due to the return of s well above sh- Only 
for the highest noise levels considered (D = 3 x 10"'^) the 
noise drives an early escape, reflected in the ISI density 
concentrated at shorter ISI lengths. 

In Fig. [24] we show the coherence measure /3 for the 
two parameter sets used for generating Class 2 and Class 
3 time series. As with the earlier analyses, obtaining fi 
for Class 1 time series is not possible due to a multi- 
peaked PSD. For Classes 2 and 3, we neglect a strong 
peak in the PSD resulting from the concatenation and 
compute /3 as in previous examples. There is a weak 
peak in /3 at intermediate noise levels, consistent with 
previous examples that exhibit MMOs related to both 
CR and SPHB. In Eqs. [I0][l2|this combination of STO 
mechanisms results from the return of s and its slower 
variation near the HB. 
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FIG. 23. Average amplitude before a spike (top) and ISI 
density (middle) for time series similar to Class 2 from the 
self-coupled FN-type {5 = 0.3) and various noise strengths. 
Bottom: ISI density for Class 3-like time series from the 
self-coupled MFN [5 = 0.17). Tsto ~ 40. 




FIG. 24. The coherence measure /3 (see Subsec. Ill Cp for the 
self-coupled FN model at two values of S. 



VII. CONCLUSION 

Attention towards oscillatory neural dynamics has ex- 
panded to include coherent subthreshold oscillatory ac- 
tivity in the voltage of some neuron classes as well as 
spikes^. Combined, the sequence of subthreshold oscil- 
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lations (STOs) followed by a spike comprises a class of 
mixed-mode oscillations MMOs. As these MMOs have 
been recognized more frequently in both experimental 
and modeling settings there is an increasing challenge 
to distinguish between different MMO-generating mecha- 
nisms. For stochastic systems there are additional routes 
for MMOs that include noise-induced oscillations and 
transients that are sustained with noise. 

In this article, we analyzed MMOs of this type both 
in a biophysical model for the dynamics of stellate cells 
and in augmented versions of phenomenological FN-type 
models. Our analysis focuses on the STO part of the sig- 
nal where the slow dynamics are prominent and noise has 
the greatest impact on the STOs. In both model types we 
observe signatures of two distinct oscillation-generating 
mechanisms: slow passage through a Hopf bifurcation 
(SPHB) and noise-induced coherent oscillations with fre- 
quencies close to that of the Hopf bifurcation. The latter 
of these appears in certain types of coherence resonance 
(CR). 

For noise levels that are in the range of what is ob- 
served experimentally, features of the underlying deter- 
ministic models can be hidden or transformed in the 
stochastic time series. A further complication in iden- 
tifying the MMO mechanism in stochastic models is that 
substantially distinct models can easily be tuned to pro- 
duce very similar time series. Then model calibration 
based on time series alone must search through a larger 
parameter space or range of stochastic models needed for 
real world data. 

We use a suite of measures to be able to distinguish be- 
tween different MMO- and STO-generating mechanisms 
in the presence of noise, differentiating between the types 
of underlying models as well as analyzing the influence 
of certain general model parameters. Identifying possible 
MMO mechanisms in this way limits the broad range of 
parameter or model space for finding appropriate MMO 
mechanisms or calibrating a biophysical model. Further- 
more this type of comparison can also identify appropri- 
ate classes of reduced models, with appropriate stochastic 
behavior over a range of parameters, to be used to ap- 
proximate a full biophysical model. Such reduced models 
are used both for simplicity and for computational speed 
within larger modeling frameworks, as, e.g., in 

The suite of these measures was chosen to focus on the 
STOs in the ISI, where noise has the greatest impact on 
the character of the MMOs. We show that the measures 
can also be used to exploit the noise to identify the route 
to the MMOs. Given the variety of routes to MMOs, 
more than one such measure is needed, and we use three 
distinct measures (interspike interval, trend of the STO 
amplitude and noise-dependent coherence) . While the fo- 
cus of these measures is on bifurcation parameters, noise 
levels, and slowly varying control parameters, we found 
that these measures also reveal information about the re- 
fractory behavior or reset in the models. Furthermore, 
we have focused on using an approach that can easily 
be applied to time series, so that it can be used in both 



experimental and simulation settings. 

We summarize the main characteristics obtained from 
these measures for different MMO generating mecha- 
nisms. 

STOs of the CR-type have distinct features in these 
measures in comparison with those dominated by SPHB: 

1. MMOs driven primarily by the deterministic be- 
havior of SPHB display a strong trend in increasing STO 
amplitude, while those driven by CR have a weaker trend 
or no trend at all, for average amplitude. 

2. For small noise ISI densities are highly concentrated 
for families driven by SPHB, in contrast to the STOs of 
the CR-type that have ISI densities with long tails. For 
larger noise values, these ISI densities are more similar 
and may depend on the return mechanism (see items 7-9 
below). 

3. The coherence measure has a clear optimal noise 
level for CR-drivcn STOs. PSDs are typically multi- 
peaked for SPHB for smaller noise, so that a coherence 
measure is not well defined for small noise. For larger 
noise values, the coherence measure for SPHB is strictly 
decreasing. 

While each model has a HB for the underlying 2D sub- 
system, differences in the criticality of these bifurcations 
can be observed depending on the underlying attracting 
deterministic dynamics. 

4. For STOs dominated by SPHB, an increase in noise 
level typically drives a greater reduction in STO ampli- 
tude before the spike for a subcritical HB than a super- 
critical HB, since the latter has attracting STOs. 

5. The underlying deterministic dynamics influence 
different behaviors of /3 as a function of noise level. Sys- 
tems with a supercritical HB may have underlying at- 
tracting behavior of small oscillations that can support 
increased coherence. This type of deterministic behavior 
is typically not seen for subcritical HB. 

6. For STOs dominated by SPHB in the subcritical 
case, even very low noise levels can drive a sampling from 
a variety of MMO families, together with a shift and 
spread of ISI values, making it difficult to distinguish 
from MMOs driven by other mechanisms. 

In addition to identifying features related to noise level 
and bifurcation or slowly varying control parameters, 
we identified a number of additional behaviors that are 
related to the reset or refractory dynamics following a 
spike. 

7. Reset and rate of increase of the control parameter 
in the simple IF model can be varied to capture the am- 
plitude and ISI behavior of the different classes of time 
series in the physiological model. However, a simple IF 
model may not have the flexibility to capture different 
behaviors of the coherence measure. 

8. For the dynamic voltage-dependent control parame- 
ters, an early escape to spiking can translate into a lower 
return value. In that case the ISI density does not shift 
with larger noise but is just less concentrated. This be- 
havior distinguishes it from IF-type models with a fixed 
reset or models where the return mechanisms are inde- 
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pendent of the spike. 

9. For return mechanisms with hmited damping of 
STOs following the spike, the coherence measure shows 
only a decrease in the coherence measure or no well- 
defined coherence over the ISI. 

One important element in these comparisons is that 
characteristics of the time series, captured by the suite 
of measures, can change in distinct ways when the noise 
level is varied. A direct way of varying noise levels in 
neuroscience experiments is by injecting a noisy current. 
Recognizing that the noise has its greatest impact in the 
ISI where the slow dynamics are prominent leads to the 
proposal that tunable extrinsic noise can be introduced 
through an applied current to mimic the effects of the 
intrinsic noise. We show that this is indeed the case with 
an appropriate scaling of the noise in the physiological 
model. 

The distinctions between model features highlighted 
by these measures have been observed in an additional 
measure related to spike clusters, repeated spikes with- 
out STOs in the ISI^^. There it was shown that the 
dynamic return mechanism of NLMFN, where an earlier 
spike can result in a reset value farther from the HB, can 
result in more robust STOs in the ISI, consistent with 
the ISI density behavior of NLMFN. The spike cluster 
frequency increases much more dramatically with noise 
for systems with SPHB-driven STOs, as compared with 
the case where STOs are CR-driven, consistent with the 
amplitude and ISI density results observed here. Also, 
small perturbations in the reset value can translate into 
a large increase in spike cluster frequency in stochastic 
systems, particularly where the STOs are SPHB-driven 
for an underlying subcritical HB. 

Our analysis in this paper is focused on a specific type 
of MMOs, namely a combination of small amplitude os- 
cillations and spikes relevant for neural systems. We have 
proposed measures that are focused on characteristics of 
MMOs that are particularly sensitive to the noise, due 
to the presence of multiple scales. This suggests that 
understanding the presence of multiple time scales and 
noise-sensitive characteristics of the underlying bifurca- 
tion structures would provide a solid basis for identifying 
measures that characterize MMOs in broader or more 
generic settings. 
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Appendix A: Abbreviations 



STO subthreshold osciUation 

MMO mixed mode osciUation 

ISI interspike interval 

PSD power spectral distribution 

IF integrate and fire 

SC stellate cell 

FN FitzHugh-Nagumo 

MFN modified FitzHugh-Nagumo 

LMFN linearly augmented modified FN 

NLMFN nonlinearly augmented modified FN 

vdP van der Pol 

HB Andronov-Hopf bifurcation 

CR coherence resonance 

SPHB slow passage through a Hopf bifurcation 



Appendix B: Equations for the seven dimensional 
biophysical model for the stellate cells 



The full seven dimensional system for the stellate cells 
(7DSC) as presented in Ref. [l! consists of one differential 
equation for the transmembrane voltage and six for six 
gating variables. Throughout this article, we omit the 
units both in the equations and the parameters. Voltage 
and reversal potentials are in mV, time in ms, all gating 
variables are unitless, the membrane capacitance C is in 



and the conductances in 



mS 



4" 
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The parameters used throughout this article are as fol- 
lows: 

C = 1; Gh ~ 1.5; Gp = 0.5; Gl = 0.5; Gk = 11; = 52; 
Eh = -20; E^^ = 55; E^ = -65; E^ = -90. (B8) 

Equivalent to Ref.l30iand what we did in the 3DSC model 
in Eq. [U we augment Eq. IB5I by the additive noise term 
^/2DT]{t). 

The following alternative equation for Ts was intro- 
duced in Ref. ^3 



dt 



1 



/ 



(H-exp(™^))^ 
5.6 



exp 



V 14 ) 



exp 



{ y+260 ' 
V 43 , 



+ 1 



(B9) 



The right hand sides of Eqs. IB7I (Eq. [51) and lBQI are very 
similar within the relevant parameter ranges of and V 
with the only noticeable deviation at very small values 
of V (between —80 and —70). The qualitative features 
of the reduced 3DSC model therefore are not altered. 
Here, we use Eq. IB7I throughout the paper, except for 
Subsec flVBi where we replace Eq. IB7l bv Eq. IBQI for the 
sake of comparability with the analysis in Ref. [S^- We 
refer to the respective footnote in Ref. [13 for a discussion. 
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Appendix C: Details for tlie measures used to characterize 
MMOs 

1. Amplitude trend before a spike 

First, the data (time series) was low-pass filtered by 
convolving it with a triangle function (e.g., Ref. Is^ ) of 
a length that is on the order of half a STO period. A 
simple algorithm then finds the spikes, chooses a window 
before each spike, removes the average and finds the local 
maxima in each window going backwards in time. The 
amplitude at each consecutive maximum before a spike 
is averaged over all ISIs in the time series. In this anal- 
ysis the choice of the low pass filter ('triangle length') is 
crucial, to avoid reduction of amplitudes or miscounting 
of the order of the maxima resulting from the time scale 
of the filter being too short or long. 



2. Power spectral distribution 

For the spiking models, spikes were removed from the 
time series by deleting data points corresponding to a 
typical spike duration (including recovery) after the re- 
spective variable {V or u) crossed a certain value. In 
the IF-type models, a similar but much shorter stretch 
of data was removed. The remaining data points were 
concatenated. Typically, PSDs computed from the full 
time series (including spikes and recovery) show signifi- 
cant power in broad frequency ranges that differ strongly 
between different models due to the different shapes of 
the spikes. Removing the spikes makes the power con- 
tribution of the STOs more clearly visible. To remove 
low frequency content in the PSDs, the average of the 
time series was removed. The routine spctrm from nu- 
merical recipes^ was used to compute an estimate of the 
PSD. The normalization is such that the total power in 
the PSD is equal to the mean squared amplitude of the 



time series. Our 'standard' PSD for the MFN models 
is obtained from time series sampled at At = 0.02 with 
a Bartlett window of size 4096 data points with over- 
lap2^. The 'standard' PSD for the SC models is obtained 
from time series sampled at At = 2.5 with the same win- 
dow function. For time series sampled at a different At, 
we scale the PSD accordingly. We scale the PSDs such 
that the frequencies are expressed in I/Tsto (see Sub- 
sec. HTC]). 

3. Coherence measure 13 

To compute the coherence measure /3 according to 
Eq. [HI we obtain the relevant values of the PSD by fit- 
ting a non-normalized Cauchy-Lorentz distribution to the 
peak of the PSD corresponding to STOs. 

The fit is generally good for intermediate and strong 
values of noise. For weaker noise values the PSDs are 
often multi-peaked making it difficult to obtain a mean- 
ingful coherence measure. The coherence measure /? has 
units of the PSD (squared unit of the considered variable 
when using Tsto as a time scale) and depends strongly 
on the sampling rate of the original time series and the 
windowsize of the PSD-estimator. It is therefore difficult 
to compare absolute values of /3 between time series with 
different units, scales and sampling rates. 
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